## Calculate atmospheric fraction of CO2
a0 = .217
a1 = .259
a2 = .338
a3 = .186
tau1 = 172.9
tau2 = 18.51
tau3 = 1.186

t <- 0:1000
af <- a0 + a1*exp(-t/tau1) + a2*exp(-t/tau2) + a3*exp(-t/tau3) # Atmospheric fraction

## Convert CO2 pulse into additional forcing
stock <- 410*7.8 # 410ppm times 7.8 GtCO2.
rf <- 5.35 * log((stock + 100*af)/stock) # Addition to radiative forcing

## Add forcing to RCPs
rcp <- read.csv("~/Box Sync/IAM_stochastic/rf_rcp2_6.txt",sep="\t",header=FALSE)
rcp$V3 <- rcp$V2
rcp$V3[rcp$V1>2019] <- rcp$V2[rcp$V1>2019] + rf[1:nrow(rcp[rcp$V1>2019,])]

write.table(rcp,"~/Box Sync/IAM_stochastic/rf_rcp2_6_plus.txt", row.names = FALSE, col.names = FALSE)